import pandas as pd
import numpy as np
import random
import matplotlib.pyplot as plt
import os
import time

# 定义模拟退火算法类
class SimulatedAnnealing:
    def __init__(self, initial_params, param_ranges, max_iterations, initial_temp, cooling_rate, data):
        self.k, self.a, self.b, self.c = initial_params  # 初始参数
        self.param_ranges = param_ranges  # 参数范围
        self.max_iterations = max_iterations  # 最大迭代次数
        self.temp = initial_temp  # 初始温度
        self.cooling_rate = cooling_rate  # 降温速率
        self.data = data  # 数据
        self.best_solution = (self.k, self.a, self.b, self.c)  # 最优解
        self.best_error, self.best_mae, self.best_rmse = self.objective_function(self.k, self.a, self.b, self.c)  # 最优误差
        self.error_history = []  # 记录每次迭代的误差

    # 目标函数
    def objective_function(self, k, a, b, c):
        P_actual = self.data['P']
        f = self.data['f']
        B = self.data['B']
        t = self.data['t']
        P_predicted = k * (f ** a) * (B ** b) * (1 + c * t)

        # 相对误差
        relative_error = np.mean(((P_actual - P_predicted) / P_actual) ** 2)
        # 平均绝对误差 (MAE)
        mae = np.mean(np.abs(P_actual - P_predicted))
        # 均方根误差 (RMSE)
        rmse = np.sqrt(np.mean((P_actual - P_predicted) ** 2))

        return relative_error, mae, rmse

    # 执行一次模拟退火迭代
    def run(self):
        for iteration in range(self.max_iterations):
            # 生成新的解
            new_k, new_a, new_b, new_c = self.generate_new_solution(self.k, self.a, self.b, self.c)
            new_error, new_mae, new_rmse = self.objective_function(new_k, new_a, new_b, new_c)

            # 比较新解与当前解
            if new_error < self.best_error:
                self.k, self.a, self.b, self.c = new_k, new_a, new_b, new_c
                self.best_error = new_error
                self.best_mae = new_mae
                self.best_rmse = new_rmse
                self.best_solution = (new_k, new_a, new_b, new_c)
            else:
                # 按一定概率接受较差解
                if self.accept_solution(new_error):
                    self.k, self.a, self.b, self.c = new_k, new_a, new_b, new_c

            # 降低温度
            self.temp *= self.cooling_rate
            self.error_history.append(self.best_error)

        return self.best_solution, self.best_error, self.best_mae, self.best_rmse

    # 生成新解
    def generate_new_solution(self, current_k, current_a, current_b, current_c):
        perturbation_k = random.uniform(-0.1, 0.1)
        perturbation_a = random.uniform(-0.1, 0.1)
        perturbation_b = random.uniform(-0.1, 0.1)
        perturbation_c = random.uniform(-0.1, 0.1)

        new_k = current_k + perturbation_k
        new_a = current_a + perturbation_a
        new_b = current_b + perturbation_b
        new_c = current_c + perturbation_c

        # 确保新参数在指定范围内
        new_k = max(self.param_ranges[0][0], min(self.param_ranges[0][1], new_k))
        new_a = max(self.param_ranges[1][0], min(self.param_ranges[1][1], new_a))
        new_b = max(self.param_ranges[2][0], min(self.param_ranges[2][1], new_b))
        new_c = max(self.param_ranges[3][0], min(self.param_ranges[3][1], new_c))

        return new_k, new_a, new_b, new_c

    # 接受新解的概率
    def accept_solution(self, new_error):
        delta_error = new_error - self.best_error
        acceptance_probability = np.exp(-delta_error / self.temp)
        return random.random() < acceptance_probability


# 记录运行开始时间
start_time = time.time()

# 读取CSV文件路径
file_path = r'D:\2024华为杯数学建模\2024C题-数据驱动与电磁损耗\C题\附件一数据拆解\Q2\Q2磁芯材料1正弦波形.csv'
if not os.path.exists(file_path):
    raise FileNotFoundError(f"文件不存在: {file_path}")
try:
    data = pd.read_csv(file_path, skiprows=1, encoding='gbk')
except UnicodeDecodeError:
    data = pd.read_csv(file_path, skiprows=1, encoding='ISO-8859-1')

data.columns = ['t', 'f', 'P', 'Unknown', 'B']  # 重命名列

# 定义参数范围
param_ranges = [
    (0.05, 0.15),  # k 的范围
    (1.5, 1.8),  # a 的范围
    (2.5, 2.7),  # b 的范围
    (-1, 1)  # c 的范围
]

# 初始参数
initial_params = [
    random.uniform(0.05, 0.15),  # k 的初始值
    random.uniform(1.5, 1.8),  # a 的初始值
    random.uniform(2.5, 2.7),  # b 的初始值
    random.uniform(-1, 1)  # c 的初始值
]

# 保存每次迭代结果
results = []

# 创建输出目录
output_dir = r'D:\2024华为杯数学建模\2024C题-数据驱动与电磁损耗\C_code\Q2_code\Q2可视化结果\Q2修正图示\线性'
os.makedirs(output_dir, exist_ok=True)

# 运行50次模拟退火算法
for run in range(50):
    # 实例化模拟退火算法
    sa = SimulatedAnnealing(
        initial_params=initial_params,
        param_ranges=param_ranges,
        max_iterations=1000,  # 设置迭代次数
        initial_temp=1000,  # 初始温度
        cooling_rate=0.95,  # 降温速率
        data=data
    )

    # 运行模拟退火算法
    best_solution, best_error, best_mae, best_rmse = sa.run()

    # 将每次结果保存到列表
    best_k, best_a, best_b, best_c = best_solution
    results.append([run + 1, best_k, best_a, best_b, best_c, best_error, best_mae, best_rmse])

    print(f"第{run + 1}次运行结果: k = {best_k}, a = {best_a}, b = {best_b}, c = {best_c}, 最小误差 = {best_error}, MAE = {best_mae}, RMSE = {best_rmse}")

    # 可视化实际值与预测值
    f = data['f']
    B = data['B']
    t = data['t']
    P_actual = data['P']
    P_predicted = best_k * (f ** best_a) * (B ** best_b) * (1 + best_c * t)

    plt.figure(figsize=(10, 6))
    plt.scatter(P_actual, P_predicted, color='blue', label='preP')
    plt.plot([P_actual.min(), P_actual.max()], [P_actual.min(), P_actual.max()], 'r--', label='y=x')
    plt.xlabel('trueP')
    plt.ylabel('preP')
    plt.title(f'prePVStrueP - run_num {run + 1}')
    plt.legend()
    plt.grid()
    plt_path = os.path.join(output_dir, f'真实P_vs_预测P_{run + 1}.png')
    plt.savefig(plt_path)
    plt.close()

    # 绘制迭代曲线
    plt.figure(figsize=(10, 6))
    plt.plot(range(1, len(sa.error_history) + 1), sa.error_history, label='best error')
    plt.xlabel('iteration')
    plt.ylabel('error')
    plt.title(f'error evolution {run + 1}')
    plt.legend()
    plt.grid()
    convergence_path = os.path.join(output_dir, f'迭代曲线_运行_{run + 1}.png')
    plt.savefig(convergence_path)
    plt.close()

# 记录运行结束时间
end_time = time.time()

# 计算并输出运行时间
execution_time = end_time - start_time
print(f"运行50次的总时间: {execution_time:.2f} 秒")

# 将所有运行结果保存到 CSV 文件
output_csv_path = os.path.join(output_dir, '50次拟合结果.csv')
result_df = pd.DataFrame(results, columns=['run_num', 'k', 'a', 'b', 'c', 'relative_error', 'mae', 'rmse'])
result_df.to_csv(output_csv_path, index=False, encoding='utf-8-sig')
print(f"所有50次拟合结果已保存到: {output_csv_path}")